####MEAN LEVEL MODELS#####
#simple models while controlling for age, period, and cohort
gss_meanPolviews <- lm_robust(att ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_meansWelfare <- lm_robust(att ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_meanrGender  <- lm_robust(att ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_meancMoral   <- lm_robust(att ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_meanImmig   <- lm_robust(att ~ genPew + year + ageGroup, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#controlled models with controls for white, protestants, and education
gss_meanPolviews_C <- lm_robust(att ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_meansWelfare_C <- lm_robust(att ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_meanrGender_C  <- lm_robust(att ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_meancMoral_C   <- lm_robust(att ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_meanImmig_C   <- lm_robust(att ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#run the same models with party id as the dependent variable 
gss_meanPolviews_pid <- lm_robust(partyid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_meansWelfare_pid <- lm_robust(partyid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_meanrGender_pid  <- lm_robust(partyid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_meancMoral_pid   <- lm_robust(partyid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_meanImmig_pid   <- lm_robust(partyid ~ genPew + year + ageGroup, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#run the same models with party id as the dependent variable 
gss_meanPolviews_pid_C <- lm_robust(partyid ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_meansWelfare_pid_C <- lm_robust(partyid ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_meanrGender_pid_C  <- lm_robust(partyid ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_meancMoral_pid_C   <- lm_robust(partyid ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_meanImmig_pid_C   <- lm_robust(partyid ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

####DIVERGENCE MODELS#####
#first, obtain the absolute residuals from the previous models
gss_Polviews$div  <- abs(gss_meanPolviews$fitted.values - gss_Polviews$att)
gss_Swelfare$div  <- abs(gss_meansWelfare$fitted.values - gss_Swelfare$att)
gss_Racgender$div <- abs(gss_meanrGender$fitted.values - gss_Racgender$att)
gss_Cultmoral$div <- abs(gss_meancMoral$fitted.values - gss_Cultmoral$att)
anes_Immig$div    <- abs(anes_meanImmig$fitted.values - anes_Immig$att)

#repeat for party id
gss_Polviews$div_pid  <- abs(gss_meanPolviews_pid$fitted.values - gss_Polviews$partyid)
gss_Swelfare$div_pid  <- abs(gss_meansWelfare_pid$fitted.values - gss_Swelfare$partyid)
gss_Racgender$div_pid <- abs(gss_meanrGender_pid$fitted.values - gss_Racgender$partyid)
gss_Cultmoral$div_pid <- abs(gss_meancMoral_pid$fitted.values - gss_Cultmoral$partyid)
anes_Immig$div_pid    <- abs(anes_meanImmig_pid$fitted.values - anes_Immig$partyid)

#same procedure for the models with demographic controls
gss_Polviews$div_C  <- abs(gss_meanPolviews_C$fitted.values - gss_Polviews$att)
gss_Swelfare$div_C  <- abs(gss_meansWelfare_C$fitted.values - gss_Swelfare$att)
gss_Racgender$div_C <- abs(gss_meanrGender_C$fitted.values - gss_Racgender$att)
gss_Cultmoral$div_C <- abs(gss_meancMoral_C$fitted.values - gss_Cultmoral$att)
anes_Immig$div_C    <- abs(anes_meanImmig_C$fitted.values - anes_Immig$att)

#repeat for party id
gss_Polviews$div_pid_C  <- abs(gss_meanPolviews_pid_C$fitted.values - gss_Polviews$partyid)
gss_Swelfare$div_pid_C  <- abs(gss_meansWelfare_pid_C$fitted.values - gss_Swelfare$partyid)
gss_Racgender$div_pid_C <- abs(gss_meanrGender_pid_C$fitted.values - gss_Racgender$partyid)
gss_Cultmoral$div_pid_C <- abs(gss_meancMoral_pid_C$fitted.values - gss_Cultmoral$partyid)
anes_Immig$div_pid_C    <- abs(anes_meanImmig_pid_C$fitted.values - anes_Immig$partyid)

#next, run the lm models using the absolute residuals as the output variables
gss_divPolviews <- lm_robust(div ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_divsWelfare <- lm_robust(div ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_divrGender  <- lm_robust(div ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_divcMoral   <- lm_robust(div ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_divImmig   <- lm_robust(div ~ genPew + year + ageGroup, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#again, same procedure for the demographic models
gss_divPolviews_C <- lm_robust(div_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_divsWelfare_C <- lm_robust(div_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_divrGender_C  <- lm_robust(div_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_divcMoral_C   <- lm_robust(div_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_divImmig_C   <- lm_robust(div_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#repeat for partyid
gss_divPolviews_pid <- lm_robust(div_pid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_divsWelfare_pid <- lm_robust(div_pid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_divrGender_pid  <- lm_robust(div_pid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_divcMoral_pid   <- lm_robust(div_pid ~ genPew + year + ageGroup, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_divImmig_pid   <- lm_robust(div_pid ~ genPew + year + ageGroup, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#repeat#2
gss_divPolviews_pid_C <- lm_robust(div_pid_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)
gss_divsWelfare_pid_C <- lm_robust(div_pid_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)
gss_divrGender_pid_C  <- lm_robust(div_pid_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)
gss_divcMoral_pid_C   <- lm_robust(div_pid_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
anes_divImmig_pid_C   <- lm_robust(div_pid_C ~ genPew + year + ageGroup + college + white + protcattend, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

#compute variables in main data frame
gss_Polviews <- gss_Polviews %>%
  mutate(sort_dv = (polviews - gss_meanPolviews$fitted.values) / gss_divPolviews$fitted.values,
         sort_iv = (as.numeric(partyid) - gss_meanPolviews_pid$fitted.values) / gss_divPolviews_pid$fitted.values,
         sort_dv_C = (polviews - gss_meanPolviews_C$fitted.values) / gss_divPolviews_C$fitted.values,
         sort_iv_C = (as.numeric(partyid) - gss_meanPolviews_pid_C$fitted.values) / gss_divPolviews_pid_C$fitted.values)

gss_Swelfare <- gss_Swelfare %>%
  mutate(sort_dv = (att - gss_meansWelfare$fitted.values) / gss_divsWelfare$fitted.values,
         sort_iv = (partyid - gss_meansWelfare_pid$fitted.values) / gss_divsWelfare_pid$fitted.values,
         sort_dv_C = (att - gss_meansWelfare_C$fitted.values) / gss_divsWelfare_C$fitted.values,
         sort_iv_C = (partyid - gss_meansWelfare_pid_C$fitted.values) / gss_divsWelfare_pid_C$fitted.values)

gss_Racgender <- gss_Racgender %>%
  mutate(sort_dv = (att - gss_meanrGender$fitted.values) / gss_divrGender$fitted.values,
         sort_iv = (partyid - gss_meanrGender_pid$fitted.values) / gss_divrGender_pid$fitted.values,
         sort_dv_C = (att - gss_meanrGender_C$fitted.values) / gss_divrGender_C$fitted.values,
         sort_iv_C = (partyid - gss_meanrGender_pid_C$fitted.values) / gss_divrGender_pid_C$fitted.values)

gss_Cultmoral <- gss_Cultmoral %>%
  mutate(sort_dv = (att - gss_meancMoral$fitted.values) / gss_divcMoral$fitted.values,
         sort_iv = (partyid - gss_meancMoral_pid$fitted.values) / gss_divcMoral_pid$fitted.values,
         sort_dv_C = (att - gss_meancMoral_C$fitted.values) / gss_divcMoral_C$fitted.values,
         sort_iv_C = (partyid - gss_meancMoral_pid_C$fitted.values) / gss_divcMoral_pid_C$fitted.values)

anes_Immig <- anes_Immig %>%
  mutate(sort_dv = (att - anes_meanImmig$fitted.values) / anes_divImmig$fitted.values,
         sort_iv = (partyid - anes_meanImmig_pid$fitted.values) / anes_divImmig_pid$fitted.values,
         sort_dv_C = (att - anes_meanImmig_C$fitted.values) / anes_divImmig_C$fitted.values,
         sort_iv_C = (partyid - anes_meanImmig_pid_C$fitted.values) / anes_divImmig_pid_C$fitted.values)

####SORTING MODELS#####
sort_gssPolviews <- lm_robust(sort_dv ~ sort_iv * genPew + year * sort_iv + ageGroup * sort_iv, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)

#social welfare
sort_gssSwelfare <- lm_robust(sort_dv ~ sort_iv * genPew + year * sort_iv + ageGroup * sort_iv, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)

#Race & Gender
sort_gssRacgender <- lm_robust(sort_dv ~ sort_iv * genPew + year * sort_iv + ageGroup * sort_iv, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)

#culture and morality
sort_gssCultmoral <- lm_robust(sort_dv ~ sort_iv * genPew + year * sort_iv + ageGroup * sort_iv, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)

#immigration attitudes
sort_anesImmig <- lm_robust(sort_dv ~ sort_iv * genPew + year * sort_iv + ageGroup * sort_iv, se_type = "stata", cluster = anes_Immig$year, anes_Immig)

####SORTING MODELS WITH DEMOGRAPHIC CONTROLS#####
sort_C_gssPolviews <- lm_robust(sort_dv_C ~ sort_iv_C * genPew + year * sort_iv_C + ageGroup * sort_iv_C +
                           college * sort_iv_C + white * sort_iv_C + protcattend * sort_iv_C, se_type = "stata", cluster = gss_Polviews$year, gss_Polviews)

#social welfare
sort_C_gssSwelfare <- lm_robust(sort_dv_C ~ sort_iv_C * genPew + year * sort_iv_C + ageGroup * sort_iv_C +
                           college * sort_iv_C + white * sort_iv_C + protcattend * sort_iv_C, se_type = "stata", cluster = gss_Swelfare$year, gss_Swelfare)

#race & Gender
sort_C_gssRacgender <- lm_robust(sort_dv_C ~ sort_iv_C * genPew + year * sort_iv_C + ageGroup * sort_iv_C +
                            college * sort_iv_C + white * sort_iv_C + protcattend * sort_iv_C, se_type = "stata", cluster = gss_Racgender$year, gss_Racgender)

#culture and morality
sort_C_gssCultmoral <- lm_robust(sort_dv_C ~ sort_iv_C * genPew + year * sort_iv_C + ageGroup * sort_iv_C +
                            college * sort_iv_C + white * sort_iv_C + protcattend * sort_iv_C, se_type = "stata", cluster = gss_Cultmoral$year, gss_Cultmoral)
#immigration attitudes
sort_C_anesImmig <- lm_robust(sort_dv_C ~ sort_iv_C * genPew + year * sort_iv_C + ageGroup * sort_iv_C +
                         college * sort_iv_C + white * sort_iv_C + protcattend * sort_iv_C, se_type = "stata", cluster = anes_Immig$year, anes_Immig)